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Abstract. We explore physically self-consistent models of dusty molecular tori in Active Galactic Nuclei (AGN) with the 
goal of interpreting VLTI observations and fitting high resolution mid-IR spectral energy distributions (SEDs). The input dust 
distribution is analytically calculated by assuming hydrostatic equilibrium between pressure forces - due to the turbulent motion 
of the gas clouds - and gravitational and centrifugal forces as a result of the contribution of the nuclear stellar distribution and 
the central black hole. For a fully three-dimensional treatment of the radiative transfer problem through the tori we employ the 
Monte Carlo code MC3D. We find that in homogeneous dust distributions the observed mid-infrared emission is dominated by 
the inner funnel of the torus, even when observing along the equatorial plane. Therefore, the stratification of the distribution 
of dust grains - both in terms of size and composition - cannot be neglected. In the current study we only include the effect 
of different sublimation radii which significantly alters the SED in comparison to models that assume an average dust grain 
property with a common sublimation radius, and suppresses the silicate emission feature at 9.7 urn. In this way we are able to fit 
the mean SED of both type I and type II AGN very well. Our fit of special objects for which high angular resolution observations 
(< 0.3") are available indicates that the hottest dust in NGC 1068 reaches the sublimation temperature while the maximum dust 
temperature in the low-luminosity AGN Circinus falls short of 1000 K. 
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1. Introduction 

Active galaxies are among the most energetic objects in the uni- 
verse. With current knowledge, these objects - although con- 
taining up to 60 different variations (Ward 2003) - can be de- 
scribed by a small number of models. The differences between 
them can be reduced to geometrical effects. The most com- 
mon model is called the Unified Scheme (Antonucci 1993) of 
Active Galactic Nuclei (AGN). According to this scheme, sev- 
eral components can be distinguished: A supermassive black 
hole (10 6 - 10 10 Mq) in the centre surrounded by an accretion 
disc, reaching from the marginally stable orbit up to several 
thousands of Schwarzschild radii. By turbulent processes, ma- 
terial is accreted and the gas in the disc is heated to several hun- 
dred thousand Kelvin. Hence, the spectrum of the emitted ra- 
diation peaks in the UV/optical wavelength range. In the outer 
part, the temperature of the disc drops to about 1000K, at the 
start of the larger and geometrically thicker dust reservoir, pro- 
ducing a characteristic bump in the IR spectral range (hereafter 
the IR bump). 
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The Unified Scheme differentiates between two different 
types of active galaxies: The spectral energy distributions 
(SEDs) of type I objects show the so-called blue bump, aris- 
ing from the direct radiation of the accretion disc mentioned 
above. Broad emission lines overlay the continuum spectral 
energy distribution. They come from the so-called Broad Line 
Region (BLR) near the centre of the object. Within this area, 
gas clouds are moving at high velocities due to the deep po- 
tential well of the central black hole and, therefore, emit broad 
spectral lines. Type II objects do not show a big blue bump and 
are characterised by narrow spectral lines, arising from orbit- 
ing gas clouds further away from the gravitational centre. After 
broad spectral lines were observed in the polarised light of a 
type II object, it was proposed that both AGN types belong to 
the same class of object with a dust reservoir located in a torus 
which blocks the direct view to the centre and the BLR within 
the opening of the torus for the case of type II objects. Polarised 
light is then produced by electrons and material above the open- 
ing angle of the dust torus. Therefore, the type of object simply 
depends on the viewing angle towards the torus. 

More direct evidence for the existence of an obscuring torus 
comes from recent observations of the prototype Seyfertll 
galaxy NGC 1068 in the mid-infrared by the VLTI with MIDI 
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(Leinert et al. 2003). Based on these interferometric data, Jaffe 
et al. (2004) distinguished two different dust components: a hot 
component in the centre with a diameter of less than 1 pc and 
a temperature of more than 800 K and warm dust within an 
elongated structure perpendicular to the jet axis with a temper- 
ature of approximately 320 K and a size of 3.4 pc times 2.1 pc 
perpendicular and parallel to the torus axis, respectively. This 
means that geometrically thick tori are consistent with the data. 

The first more detailed radiative transfer simulations for the 
case of dusty tori were carried out by Pier & Krolik (1992). 
They used the most simple dust configuration to describe the 
observations available at that time, consisting only of spec- 
tral energy distributions - due to the large distances (several 
tens of Mpc) and small sizes (less than lOOpc) of these ob- 
jects: a cylindrical shape with a cylindrical hole in the middle. 
Although pointing out that dust - in order to avoid destruction 
by hot gas - must be contained in clumps, they applied a ho- 
mogeneous dust distribution. Characteristics of their modelling 
are the small sizes of a few pc in diameter and very high dust 
densities. The small amount of cold dust leads to a too narrow 
dust temperature range, resulting in spectra insufficiently ex- 
tended towards longer wavelengths, compared to observations. 
Granato & Danese (1994) preferred larger wedges (tens to hun- 
dreds of pc in diameter) and thereby solved the problem of the 
too narrow IR bumps compared to observations. Further dif- 
ficulties when comparing calculations with data emerge from 
the so-called silicate feature problem. The silicate feature is a 
resonance in the spectral energy distributions at 9.7 um, arising 
from the stretching modes in silicate tetrahedra. For optically 
thin configurations, this feature is seen in emission. Looking 
at optically thick tori, two cases have to be distinguished: The 
feature appears in emission, if the temperature decreases along 
the line of sight away from the observer and an absorption fea- 
ture is seen for the case of increasing temperature along the 
line of sight. Between those two cases, at an optical depth 
i~9.7f.tm ~ 1, a mixture of both - so-called self-absorption - 
is visible (Henning 1983). Therefore, for typical torus geome- 
tries, the feature is expected to arise in absorption for type II 
objects and in emission in type I sources. Up to now, no ob- 
servation of a type I object has shown the feature in emission. 
In both studies mentioned, a significant reduction of the sil- 
icate emission-feature could be found for the case of type I 
galaxies. However, a lot of fine tuning of the parameters was 
needed. Granato & Danese (1994) point out that small sili- 
cate grains are selectively destroyed by radiation pressure in- 
duced shocks in the inner part of the torus close to the cen- 
tral source. Therefore, they set up a depletion radius for small 
silicate grains and are able to reduce the silicate emission at 
9.7mn. Manske et al. (1998) succeeded in avoiding the feature 
by using a toroidal density distribution comparable to the one 
of Granato & Danese (1994) with a large optical depth in com- 
bination with a strong anisotropic radiation source, which is a 
more physical approach to describe the accretion disc. Another 
more realistic way to cope with the silicate feature problem was 
proposed by Nenkova et al. (2002). They introduce a clumpy 
structure again within a flared disc geometry. Their radiative 
transfer simulations show that within a large range of values 
of their parameters, a reduction of the emission feature is fea- 



sible. Another possibility to investigate AGN torus models is 
the simulation of polarisation maps and their comparison with 
observations. For a first approach see Wolf & Henning (1999). 

In this first paper of a series, we present radiative trans- 
fer calculations for a hydrostatic model of dusty tori. We use 
the so-called Turbulent Torus Model, introduced by Camenzind 
(1995). The great advantage of this model arises from the fact 
that the dust density distribution and the geometrical shape of 
the torus do not have to be set up arbitrarily and independent of 
each other, but both result from physically reasonable assump- 
tions. Our calculations combine structure modelling and fully 
three-dimensional Monte Carlo radiative transfer simulations. 
The following chapter will summarise the structure modelling 
together with all ingredients needed to perform the radiative 
transfer simulations. 

In a further step towards more physical models of AGN 
dust tori, the homogeneous dust distribution will be replaced by 
a clumpy dust model. Later, we will start with hydrodynamic 
simulations in combination with radiative transfer calculations, 
in order to gain information about the temporal evolution of 
spectral energy distributions and surface brightness maps of 
dusty tori. 

Chapter 3 deals with the results we obtained with the first 
approach. Spectral energy distributions and surface brightness 
maps are calculated and discussed for part of the simulation 
series we conducted. In chapter 4 we discuss comparisons of 
our model results with large aperture data as well as with very 
high spatial resolution data. 

2. The model 

2.1. The Turbulent Torus Model 

It is known that many spiral galaxies harbor young nuclear star 
clusters (e.g. Gallimore & Matthews 2003; Walcher et al. 2004; 
Boker et al. 2002, 2004). Especially during late stages of stel- 
lar evolution, stars release large amounts of gas through stellar 
winds and the ejection of planetary nebulae. With increasing 
distance from the star, dust forms from the gas phase. Due to 
the origin of the dust from ejection processes of single stars, 
we expect the dust torus to consist of a cloudy structure. As 
even current interferometric instruments are not able to resolve 
single clouds of the dust distribution, they are assumed to be 
small. For the sake of simplicity, we treat the torus as a contin- 
uous medium in this paper. The dust clouds are ejected by the 
stars and therefore take over the stellar velocity dispersion cr t . 
As we neglect interaction effects between these clouds, we can 
assume that they possess the same velocity as the stars, called 
the turbulent velocity v t (v t « cr*). 

We are dealing with stationary models only. Therefore, the 
presence of continuous energy feedback by supernovae explo- 
sions as well as continuously injected mass by stellar winds is 
not taken into account. The dust clouds are embedded into hot 
gas, which leads to the destruction of dust grains by sputtering 
effects. Therefore, dust has a finite lifetime and it is likely that 
the composition differs from the one found in interstellar dust 
within our galaxy (see Sect. 2.3). The effects of magnetic fields 
on charged dust particles are neglected. The radiation pressure 
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of the extreme radiation field of the central AGN exerted on the 
dust grains only affects the very innermost part of the dust den- 
sity distribution and is therefore not taken into account within 
our simulations. 

Due to their environment, the dusty clouds are subject to 
gravitational forces. Together with centrifugal forces - because 
of the rotation of the central star cluster - an effective poten- 
tial can be formed. In the following, the three main compo- 
nents will be discussed: First, the gravitational potential of the 
supermassive black hole. As the torus extends to several hun- 
dred thousands of Schwarzschild radii from the centre, we can 
treat the black hole as a point-like source of gravitation with a 
Newtonian potential, given in cylindrical coordinates: 

GM BH 



, x 10' 



4>bh(R,z) = - 



(1) 



The second contribution we take into account is the gravita- 
tional potential of the nuclear star cluster, dominating the con- 
tribution by stars in the galaxy core. From high resolution ob- 
servations of the centres of galactic nuclei (e. g. Gebhardt et al. 
1996) a surface brightness distribution profile can be derived. 
This can be well fitted by a two-component power law, the 
so-called Nuker law, given by the following equation (Lauer 
et al. 1995; Byun et al. 1996; Gebhardt et al. 1996; Faber et al. 
1997a): 



I(R,z) = 2^I C 



Rc 



y/R 2 + z 2 



1 + 
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where R c is the core radius (or break radius), ft the slope of the 
surface brightness distribution outside R c , 6 the slope inside R c , 
a characterises the width of the transition region and I c is the 
surface brightness at the location of R c . 

When assuming a constant mass-to-light ratio, this distribu- 
tion mirrors the projected mass density distribution directly. By 
applying a least-square fitting procedure, one can obtain the 3D 
density profile, for which an analytical expression can be given. 
This is the so-called Zhao density profile family (Zhao 1997): 



4nR 3 r 
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(3) 



where the parameters a, fi and y are complicated functions of 
the parameters a, j3 and y of the Nuker law. The quantity C is a 
constant, depending on these parameters. 

By solving the Poisson equation, the gravitational potential 
of the nuclear star cluster can be obtained. Throughout this pa- 
per, all simulations use the special case of the Hernquist profile 
(Hernquist 1990). According to observed profiles of cores of 
spiral galaxies, the exponent of the density distribution within 
the core radius is close to 1 in most cases (corresponding to the 
Hernquist profile). Recent simulations of dark matter halos by 
Williams et al. (2004) - which belong to the same family of 
stellar systems - also give an upper limit for the exponent of 
the density distribution within the core radius between 1 .5 and 



0.5 




R[pc] 



Fig. 1. Comparison of the contributions of the three components of 
the effective potential for the parameters of our standard model in the 
equatorial plane. The solid line corresponds to the effective potential, 
the dotted line to the BH-potential, the dashed line to the potential of 
the nuclear star cluster and the dashed-dotted line to the centrifugal 
potential. 



1 .7 in the inner part. 

With the parameters a = 1, ft = 4 und 6 - 
needed gravitational potential: 
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The third part of the effective potential is the centrifugal po- 
tential. Apart from the turbulent velocity, the stars - and there- 
fore the dust clouds - possess an additional velocity component 
from the rotational motion around the gravitational centre. Due 
to the lack of observational data, we adopt the following distri- 
bution of specific angular momentum for such a stellar system: 

I R V 

./'spec = VG(M BH + M t (R T ))Rj — . 

Here, /? T is the torus radius, defined as the location, where the 
effective potential reaches its minimal value. This corresponds 
to an equilibrium of centrifugal and gravitational forces in our 
solution, where the angular momentum attains the Keplerian 
value. Therefore, a distribution of the specific angular momen- 
tum is needed, which fixes the value at the location of the torus 
radius to the Keplerian value. We used a less steep decline of 
the specific angular momentum beyond the torus radius than 
for the Keplerian case, which is in agreement with expecta- 
tions for elliptical star clusters. Such toroidal gas distributions 
are expected to be unstable on time scales larger than the dy- 
namical time at the torus radius. The torus is however steadily 
refilled so that a quasi-stationary configuration results. 
By integrating over the centrifugal force - under the assump- 
tion of constant angular velocity on cylinders - we obtain the 
centrifugal potential: 

MR,z)-- G(M ™: MM) (Ar 1 '- (6) 
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The effective potential and its three single components are 
shown in Fig. 1 for the case of our standard model (described 
in Sect. 2.2) in the equatorial plane. As can be seen, the 
effective potential (given by the solid line) in the outer part 
of the torus is dominated by the potential of the nuclear 
star cluster (dashed line). Moving closer to the centre, the 
black hole potential (dotted line) wins, which leads to a 
further bending down of the potential, before the centrifugal 
potential (dashed-dotted line) dominates. The latter is the 
only repulsive term and leads to a centrifugal barrier close to 
the centre, preventing the models to be gravitationally unstable. 

The force given by the effective potential has to be bal- 
anced by the turbulent pressure in hydrostatic equilibrium. This 
pressure results from the turbulent motion of the dust clouds 
and can be expressed by the following (isothermal) equation of 
state: 



P(pd) = Pd v t 



(7) 



After inserting this into the equation for the hydrostatic equilib- 
rium and solving the differential equation under the assumption 
that v t is constant, we obtain the dust density distribution: 



Pd = P A exp 



GM t M BH /M« 
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V fR c | ^ 
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Here, p° is a measure for the total dust mass. From this equa- 
tion, some basic features of the model can be derived: First 
of all, equipotential surfaces have the same shape as surfaces 
of constant density, so-called isopynic surfaces. The shape of 
these surfaces only depends on the reduced parameter Mbh/M* 
and the reduced coordinates R/R c , R/Rj and z/R c and the ex- 
ponent of the angular momentum distribution y. As mentioned 
above, the torus radius is connected with the minimum of the 
potential, the maximum of the density and the maximum of 
pressure and lies within the equatorial plane, due to symmetry 
considerations. 

2.2. Parameters of the model 

The introduced standard model will be used to represent the 
core of a typical Seyfert galaxy. Therefore, we calculated the 
mean black hole mass and mean luminosity of a large sam- 
ple of Seyfert galaxies, given in Woo & Urry (2002). We find 
mean values of <M B h) = 6.6- 10 7 M o and (Ldisc) = 1.2- 10 n L Q . 
After calculating a set of simulations with different dust den- 
sity distributions, leading to different density contrasts between 
the maximum density and the outer part of the distribution, we 
found that a relatively shallow density distribution is needed 
in order to obtain SEDs extending to wavelengths longer than 
10-20 Lim. To achieve this, a relatively steep radial distribution 
of specific angular momentum is required (y = 0.5 leads to the 
steepest - but still stable - distribution). But this gives only a 
fairly small effect. It is more efficient to increase the core ra- 
dius, which leads to a broadening of the density distribution, 



Table 1. 

model. 



Parameters used for the simulation of our mean Seyfert 



Parameter 


Value 




6.64 • 10 7 M Q 




75 pc 


M, 


2.0 • 10" M Q 


-^disc 


1.23- 10" L Q 


Aiisc/^Edd 


6% 


Ri 


5pc 


M dust 


5.79 • 10 5 M e 


T9.7|im 


2.0 


^out 


225 pc 


Vt 


201 km s 1 


r 


0.5 


d 


45Mpc 



due to the flattening of the potential of the star cluster in the 
outer part of the torus. A change of the total mass of the stars 
mainly scales the potential of the stellar cluster and therefore 
finally the density distribution. By taking this into account, we 
chose the core radius for our mean model to be 75 pc with a 
mass of the nuclear star cluster of 2 • 10 9 M . 

To constrain the other parameters, one can use the fact that 
not all of them are independent from each other. Comparable to 
the fundamental plane of global galaxy parameters, the cores of 
ellipticals and bulges of spiral galaxies span the so-called core 
fundamental plane, relating the parameters core radius R c , the 
central surface brightness fio and the central velocity dispersion 
<x* (Kormendy 1987; Faber et al. 1997b). These studies give 
several relations between relevant parameters. The relation be- 
tween the core radius and the central velocity dispersion of the 
stars is interesting for us: 



201 



75 pc 



km 

s 



(9) 



As already mentioned above, the dusty clouds are produced 
by the stars and take over their velocity dispersion. Therefore, 
we choose the turbulent velocity of the clouds v t to be equal 
to the velocity dispersion of the stars in the central region of 
the galaxy. According to equation (9) the value v t = 201 kms -1 
is obtained. Interactions between the different phases (dust, 
molecular and hot gas) within the torus and between other dust 
clouds, which may alter their velocity distribution are not taken 
into account at present. The outer radius of the torus is cho- 
sen to be three times the core radius of the nuclear star cluster, 
where the density distribution has already reached an approxi- 
mately constant value. An upper limit for the total mass of the 
nuclear star cluster is determined by the Jeans equation, which 
relates structural with dynamical parameters for an isothermal 
model of the central stellar distribution under the assumption of 
a constant velocity dispersion. The remaining two parameters 
Rj and y are chosen such that the dust density distribution re- 
mains flat enough to obtain IR bumps, which are broad enough 
to compare with data and to obtain a torus model with an inner 
radius small enough to allow the majority of the grain species 
to reach their sublimation temperatures. Finally, the mass of the 
dust, enclosed into the torus, is determined by the depth of the 
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Fig. 2. Visualisation of our standard model: a) Isopynic line plot (logarithmically equidistant) in a meridional plane of the torus. The numbers 
give the dust density in kg irr 3 of the corresponding isopynic line, b) Dust density distribution in a meridional plane. 



silicate feature at 9.7 ^m. We chose Md us t = 5.79 • 10 5 Mq (cor- 
responding to an optical depth of T9.7 = 2.0), which yields 
depths comparable to observed spectra. 

The parameters we assumed for our standard model are 
summarised in Table 1 and the resulting model is shown in 
Fig. 2. As can be seen, the torus has a relatively dense core 
and the density distribution gets shallow in the outer part. The 
maximum is reached in the equatorial plane (due to symmetry 
considerations) at the torus radius Rj and the models possess a 
more or less pronounced cusp at the inner end of the torus. The 
dust-free cone with quite steep walls will be called funnel in 
the following. The turbulent velocity always tends to increase 
the height of the torus, while the angular momentum tends to 
flatten it. 

2.3. Dust properties 

Although the very extreme physical conditions near the central 
source make it likely that the dust composition may be altered 
compared to interstellar dust in our own galaxy, we assume - 
for the sake of simplicity and comparability - a typical dust ex- 
tinction curve of interstellar dust. We use the model of Mathis 
et al. (1977) which fits well the observed extinction curve of 
the interstellar medium. 

They found a number density distribution (MRN- 
distribution hereafter), which is proportional to the grain ra- 
dius to the power of - 3.5 and a range of grain radii between 
0.005 mn and 0.25 \im. We use 15 different grains 1 , containing 



a mixture of 5 silicate and 10 graphite grains - each species 
with different grain radii according to the MRN grain size dis- 
tribution - and a mass fraction of 62.5% silicate and 37.5% 
graphite. Because of the anisotropic behaviour of graphite, two 
different sets of optical constants are necessary: one where the 
electric field vector oscillates parallel to the crystal axis of the 
grain (5 grains) and another one, where it oscillates perpen- 
dicular to it (5 grains). For them, the standard 3-3 approx- 
imation is used, which is valid exactly in the Rayleigh limit 
and reproduces extinction curves reasonably well, as shown 
by Draine & Malhotra (1993). Optical data (refraction indices) 
are taken from Draine & Lee (1984), Laor & Draine (1993) 
and Weingartner & Draine (2001) and the scattering matrix el- 
ements and coefficients are calculated under the assumption 
of spherical grains using the routines of Bohren & Huffman 
(1983) within the radiative transfer code MC3D, which will be 
discussed in Sect. 2.5. 

In addition, a comparison was made between a split into 
15 different grains and the usage of only one grain, but with 
weighted mean values for the extinction, absorption and scat- 
tering cross sections and the Stokes parameters, which describe 
the intensity and the polarisation state of the photon package. 
The weights result from the abundance and the size distribu- 
tion of the respective material (for details about the averaging 
see Wolf (2003a)). The resulting extinction curves of these two 
cases are shown in Fig. 3a, where the solid line corresponds 
to the splitted model and the dashed line to the model with 
mean characteristics. For every dust species (type of material 
and grain size) the temperature distribution is determined sep- 



1 The number of dust grains increases the computational time for 
simulating the temperature distribution, therefore we tried to minimise 
the number of grains. Simulations, where we doubled the number of 



grains showed no significant differences in the spectral energy distri- 
butions. 
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Fig. 3. Input parameters of our simulations: a) Course of the mass extinction coefficient plotted against the wavelength for the averaged dust 
model (dashed line) and single grains (solid line): composition of silicate and graphite grains with varying grain radii (for details see text). The 
dotted straight line should clarify the offset in the continuum extinction curve, b) SED of the central accretion disc, used as the primary energy 
source. 



arately. In Fig. 3a one can clearly see the two silicate features 
at 9.7 urn and 18.5 um. The extinction curves drop down quite 
steeply at longer wavelengths and an offset can be recognised 
around the wavelengths of the two features between the contin- 
uum spectra shortward of 10 um and long ward of 20 urn, em- 
phasised by the dotted line fragment. 

2.4. Primary source 

In our modelling, the dust in the torus is solely heated by an ac- 
cretion disc in the centre of the model space. In most of the sim- 
ulations we use a point-like, isotropically radiating source with 
a luminosity of 1.2 • 10 11 L Q in our standard model, which cor- 
responds to roughly 6% of the Eddington luminosity. The spec- 
tral energy distribution (SED) is composed of different power 
laws and a Planck curve as follows. In the ultraviolet to op- 
tical wavelength regime, we use the composite mean quasar 
spectrum from Manners (2002). For wavelengths longer than 
Lyman alpha, it includes data from 2200 quasars of the Sloan 
Digital Sky Survey (Vanden Berk et al. 2001), which results in 
a power law with a spectral index of a v = -0.46 1 . Shortwards 
of Lyman alpha, data from radio-quiet quasars (Zheng et al. 
1997), taken with the HST Faint Object Spectrograph, is used. 
The data leads to a spectral index of a v = - 1 .8. For wavelengths 
longer than 10 (im, we use a decline of the spectrum according 
to the Rayleigh-Jeans branch of a Planck curve with an effec- 
tive temperature of 1000 K, reflecting the smallest temperatures 
expected from the thermal emission of the accretion disc. At 
shorter wavelengths, the spectrum decays with a spectral index 
of a v = -3. Due to the lack of observational data in this wave- 
length range, the index is chosen in agreement with theoretical 
modelling of accretion disc spectra performed by Hubeny et al. 



2 a v is defined as F v oc V" 



(2000). Taking this shape and normalising to the chosen bolo- 
metric luminosity of the accretion disc yields the input spec- 
trum shown in Fig. 3b. 

2.5. The radiative transfer code MC3D 

MC3D is a three-dimensional continuum radiative transfer 
code based on the Monte Carlo approach (Fischer et al. 
1994; Wolf et al. 1999). It is able to manage arbitrary three- 
dimensional dust and electron configurations and different pri- 
mary sources of radiation. Diverse shapes and species of dust 
grains can be implemented as well. Working in the particle pic- 
ture, the bolometric luminosity of the central source is divided 
into photon packages, so-called weighted photons. They are 
fully described by their wavelength A and their Stokes vector, 
describing the intensity and the polarisation state of the parti- 
cle. According to a provided radiation characteristic and spec- 
tral energy distribution of the source, such a weighted photon 
is emitted. Within the dust or electron configuration, absorption 
and scattering take place, where Thomson scattering, Rayleigh 
and Mie scattering is applied to the weighted photon. In a last 
step, these photons are detected after leaving the convex model 
space. 

According to this procedure, temperature distributions as 
well as - in a second step - spectral energy distributions and 
surface brightness distributions can be calculated. Temperature 
distributions have to be calculated using real radiative trans- 
fer. In order to minimise CPU time and get smoother surface 
brightness distributions, a raytracer (also inherent in MC3D) 
can be used to obtain SEDs and maps. 

Details about MC3D and the built-in features can be found 
in Wolf et al. (1999), Wolf & Henning (2000), Wolf (2001) and 
Wolf (2003b). 
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Fig. 4. Comparison of MC3D (solid line) with the code of Granato & 
Danese (1994) (dashed line). 

Pascucci et al. (2004) tested MC3D for 2D structures 
against various other continuum radiative transfer codes, where 
they found good agreement. We also performed a direct com- 
parison for the special case of AGN dust tori with the simula- 
tions of Granato & Danese (1994), one of the standard torus 
models to compare to, calculated with his grid based code. The 
results are displayed in Fig. 4, where the solid line corresponds 
to the calculation with MC3D and the dashed line to the origi- 
nal spectrum of Granato. The underlying model is a wedge ge- 
ometry with dust of constant density. The ratio between outer 
and inner radius is set to 100 with a half opening angle of the 
dust-free cone of 40 ° and an optical depth of 0.3 at 9.7 um (for 
details of the model setup see Granato & Danese (1994)). The 
original spectrum of Granato's code is taken from the library of 
SEDs described in Galliano et al. (2003). 

As can be seen in Fig. 4, good agreement between the two 
codes is found with deviations smaller than 30% within the rel- 
evant wavelength range. The reason for the remaining differ- 
ences between the two SEDs is due to the fact that only stan- 
dard calibrations of the code-parameters were used and no fine 
tuning was performed. Minor differences of the model setup are 
responsible for the remaining discrepancies among the SEDs. 

3. Results and discussion 

First, the use of single dust grains is motivated, before we in- 
troduce our model by presenting temperature distributions and 
an inclination study. After that, several parameters of our stan- 
dard model have been varied and the resulting effects will be 
discussed. 

If not noted otherwise, all of the SEDs shown are pure 
dust re-emission spectra. The direct radiation from the central 



source is always omitted, as we are only interested in the mid- 
infrared wavelength range, where the torus emission dominates 
the SED. Optical depths are given for lines of sight within the 
equatorial plane from the centre outwards. 

3.1. Averaged dust mixture versus single grains 

Detailed Monte Carlo simulations are always very time con- 
suming. In order to reduce the CPU time, an averaged dust 
grain mixture (Wolf 2003a) is often considered instead of a real 
mixture of single grains. However, temperature distributions 
have to be computed for every single grain. Therefore, com- 
putation time strongly scales with the number of dust grains 
considered. As Wolf (2003a) showed for the case of proto- 
planetary discs and we also confirmed for AGN tori, an av- 
eraged dust model works perfectly well, if the inner radius of 
the disc/torus is larger than the sublimation radii of the grains 
(see discussion later). But for most of our models, this is not 
the case, as the inner radii of the tori are often determined by 
sublimation themselves. These sublimation radii follow an ap- 
proximative formula for the local thermodynamic equilibrium. 
It only takes emission of the central source and extinction by 
dust species with smaller inner radii into account. Although the 
resulting maximal temperatures for the single grains differ by 
up to 10% relative to the assumed sublimation temperatures of 
silicate and graphite grains, this seems to be an adequate pro- 
cedure, as sublimation temperatures are not very well known. 
With this new dust model, the situation looks much different, 
as shown in Fig. 5, although the overall shape of the extinction 
curve looks similar compared to the one of the averaged grain 
model (see Fig. 3a). The typical shape of the resulting SEDs 
will be described in more detail in Sect. 3.3. In the following, 
we address just the major differences between these two cases. 
The solid line corresponds to our standard model (see Fig. 5) 
with a dust mixture containing 15 different grains (3 different 
grain species - silicate and two orientations of graphite grains 
- with 5 varying grain radii each), while the dotted line is cal- 
culated for the same torus parameters, but the dust model is 
given by a single grain type with averaged properties concern- 
ing chemical composition and grain size. Using the averaged 
grain model, the silicate features are much more pronounced 
compared to the other case. The same happens with the off- 
set in the SED shortwards of the 9.7 um feature. The reason 
for this is the varying sublimation temperatures of the different 
dust species. Silicate grains sublimate at about 1000 K, while 
graphite grains sublimate around 1500K. Smaller particles are 
heated more efficiently than larger grains. Therefore, smaller 
grains reach the sublimation temperature earlier and the inner 
radius of their density distribution must be larger than the one 
for the distribution of large grains. This effect leads to a layer- 
ing of grain species and sizes, where silicate grains are further 
out (in the radial direction) than graphite grains and the dis- 
tribution of smaller grains is shifted outwards compared to the 
distribution of larger grains. Neglecting this layering and using 
averaged properties within one grain leads to the effects shown 
in Fig. 5. The reason for these effects is that some of the dust 
species, represented by the grain, are heated to temperatures 
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Fig. 5. SEDs at inclination angles of ° and 90 ° for the case of averaged dust properties (dotted line) and dust splitted in grains of different 
species and sizes (solid line). The horizontal line segments in the left panel indicate the offset of the continuum SEDs of the individual curves. 



higher than their sublimation temperatures. This is an unphys- 
ically treatment of these components. Especially affected are 
the smallest silicate grains, which are the most abundant ones 
in our dust model. The effects are mainly more explicitly visi- 
ble characteristics of the dust extinction curve properties of sil- 
icate grains - especially the offset in the continuum spectrum 
near the feature (demonstrated by the horizontal lines in Fig. 5) 
and the strength of the feature itself. Reaching only the permit- 
ted smaller temperatures, these effects are partially balanced by 
the contribution of the other grain types. Using the minimum of 
sublimation temperatures of the mean grain components leads 
to a Wien branch shifted to longer wavelengths. 

Being mainly interested in the mid-infrared wavelength 
range, we decided to split the dust model in single grains for 
all of the models shown in this paper. All species obey the 
same density distribution, according to physical reasoning (see 
Sect. 2.1) and the mass fraction of the individual species, de- 
rived from the power law number density distribution. 

The inner radius of the distribution of every dust species is 
calculated using the formula for the local energy equilibrium 
and taking extinction by all other dust species into account. 
Spheres with these radii around the central source are left va- 
cant concerning these particular dust species. The rest of the 
calculated equipotential surfaces are filled with dust, according 
to the dust model used. 

3.2. Temperature distribution 

Fig. 6 shows the temperature distribution of our standard 
Seyfert model for the smallest silicate grain component, which 
is also the most abundant one (see Sect. 2.3). In the left panel, 
radial temperature distributions for an equatorial line of sight 



(solid line), an inclination angle of i = 30 ° (dotted line) and 
i = 20 ° (dashed line) are plotted. With increasing inclination 
angle, the intersection of the respective line of sight with the 
torus moves to larger distances from the central source (com- 
pare to Fig. 8b). Therefore, the maxima of the temperature dis- 
tributions move further out and the maximum temperature de- 
creases as well. The smaller the inclination the smaller the dust 
density at the point of the intersection of the line of sight with 
the torus. This causes shallower temperature distributions with 
decreasing inclination. 

The right panel of Fig. 6 shows various cuts in ^-direction 
through the temperature distribution at radial distances from 
the centre of 2pc, 5 pc and 20 pc. A minimum curve is expected 
because of the direct illuminated - and therefore hotter - funnel 
walls. In the outer part of the torus, the temperature contour 
lines get spherical symmetric, as can be seen by the flattening 
of the curve at 20 pc distance from the centre. 

Fig. 7 displays the temperature distributions within the 
equatorial plane for different kinds of grains. In the left panel, 
the layering of different sized silicate grains can be read off. 
The smallest grains (given by the solid line in Fig. 7a) are 
heated more efficiently and, therefore, their distribution pos- 
sesses the largest sublimation radius. Fig. 7b compares the tem- 
perature distributions within the equatorial plane for the small- 
est grains of the three dust species. The solid line corresponds 
to silicate grains, the dotted and dashed curve to graphite with 
the two orientations of the polarisation vector relative to the 
optical axis. Silicate grains are heated more efficiently than 
graphite grains. 
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Fig. 6. Temperature distributions of the smallest silicate grain component: a) Radial temperature distributions. The solid line corresponds to 
an equatorial line of sight, the dotted line to an inclination angle of i = 30° and the dashed line to = 20°. b) Temperature distributions in 
^-direction. The solid line corresponds to a radial distance to the centre of 2 pc, the dotted line to 5 pc and the dashed line to 20 pc. 




Fig. 7. Comparison of the temperature distributions of different dust species within the equatorial plane: a) Temperature distributions for 
different sized silicate grains: a = 0.005 urn (solid line), a = 0.013 |im (dotted line), a = 0.035 |im (dashed line), a = 0.094 urn (dashed-dotted 
line), a = 0.25 urn (dashed-dotted-dotted-dotted line), b) Temperature distributions for different dust species: the smallest silicate component 
(solid line), the smallest graphite component with E || c (dotted line) and the smallest graphite component with Elc (dashed line). 



3.3. Inclination angle study for SEDs 

In Fig. 8a we show an inclination study for SEDs. The upper- 
most curve belongs to an angle of °, which corresponds to an 
extreme Seyfert I case, where we see the torus face-on. Typical 
inclination angles of Seyfert I galaxies are between 10°and 
20°. Then, the inclination angle is increased up to 90°, the 



edge-on case, in steps of 10°. Some of the inclinations are vi- 
sualised in Fig. 8b for the case of our standard model. Only dust 
re-emission SEDs are shown, which result in the so-called IR 
bump. It arises from dust at different temperatures within the 
toroidal distribution. The two characteristic silicate features at 
9.7 urn and 18.5 \im are present in the SED. The 18.5 \im fea- 
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Fig. 8. Inclination angle study for our standard model: a) SEDs of the dust re-emission. The inclination varies from 0° (upper curve) to 90° 
(lowermost curve) in steps of 10 ° . b) Illustration of some inclination angles of lines of sight frequently used for our standard model. 



ture is less pronounced and partially hidden in the global max- 
imum of the IR bump. Just looking at the continuum spectrum 
and neglecting the silicate bands for a second, an offset in the 
flux is visible. This offset (explained in Sect. 3.1 and visualised 
in Fig. 5) is again due to the offset in the extinction curve (see 
Fig. 3a) and especially visible for the case of silicate grains. 

For i = 0°,we directly see the torus and the primary source 
and, therefore, expect the silicate feature in emission. Moving 
to higher inclination angles, the extinction to occur along the 
line of sight, mainly caused by cold dust in the outer part of 
the torus, increases. Less and less of the hotter, directly illu- 
minated parts of the torus can be seen directly. Therefore, the 
silicate feature changes from emission to absorption. With ris- 
ing extinction along the line of sight, the increase of the SEDs 
at small wavelengths is shifted to longer wavelengths and the 
IR bump gets more and more narrow. The increase at small 
wavelengths can be explained by a Planck curve with the high- 
est occuring dust temperatures (at the very inner end of the 
cusp) altered by rising extinction along the line of sight. For our 
modelling, the torus gets optically thin between 40 and 50 urn 
(due to the steep decrease of the extinction curve, compare to 
Fig. 3a). Therefore, all curves coincide and the dust configura- 
tion emits radiation isotropically for longer wavelengths. 

When we look at the extinction in the visual wavelength 
range plotted against the inclination angle, A v reaches about 
5 at an inclination angle of 10°. This means that the open- 
ing angle of our torus model is too small compared to simple 
estimates of the torus opening angle using the ratio between 
Seyfert I and II galaxies, e. g. Osterbrock & Shaw (1988) find a 
half opening angle between 29°and 39°. 



3.4. Inclination angle study for surface brightness 
distributions 

In the following section, an inclination study for surface bright- 
ness distributions is presented. The underlying model is our 
standard model introduced earlier. As wavelength we chose 
A - 13.18 um, which lies in the continuum spectrum outside 
the silicate feature and still belongs to the wavelength range 
covered by MIDI (Leinert et al. 2003). Using the Wien law 
as a rough estimate of the temperature, where the dust emits 
the maximal radiation at this wavelength, we get about 230 K. 
Looking at the steep radial temperature distributions of the var- 
ious grains (see Fig. 7b), we expect the radiation to come from 
the direct vicinity of the inner radius of the torus. The result- 
ing surface brightness distributions for inclination angles of °, 
30 °, 60 ° and 90 ° are shown in the upper row of Fig. 9. The 
images are given in a linear colour scale emphasising areas 
of maximum surface brightness. The axis labeling denotes the 
distance to the centre of the AGN in pc. When we consider 
flux scales, the object is always assumed to be at a distance of 
45 Mpc, which is a typical value for Seyfert galaxies, which 
will be observed with MIDI. 

Looking at the image with an inclination angle of 90 °(the 
last image in the first row of Fig. 9), the most central part of the 
torus funnel is visible. As already mentioned above, the dust 
species sublimate at different distances from the central source. 
Within these inner radii, no such grains can survive. Therefore, 
a layering of dust grains exists in the vicinity of the cusp. This 
is partly visible here. While the very inner end of the torus ap- 
pears quite faint due to the depletion of many grain species 
(especially the most abundant small silicate grains and also the 
smallest graphite grains), a semilunar feature is visible further 
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Fig. 9. Surface brightness distributions of our standard model with an isotropically radiating central source (first row) and with a |cos 9\ - 
radiation characteristic (second row) at inclination angles of 0°, 30°, 60° and 90° given with a linear colour scale ranging from to the 
maximum value of the respective image. 



out. At this point, the dust density increases strongly, because 
the sublimation radii of small silicate grains are reached. They 
provide the highest mass fraction of all dust species (for this 
dust model, the two smallest silicate grain populations carry 
about 43% of the total mass). Apart from this, an x-shaped 
emission area appears. Here we see the directly illuminated 
funnel walls. As the dust density peaks further out at a dis- 
tance of 5 pc in the equatorial plane, the dust density next to 
the funnel (outer equipotential lines) is quite low. Hence, only 
an x-shaped feature remains from the emission of the walls, 
which arises from a summation effect along the lines tangential 
to these walls of the funnel. The two maxima are connected via 
an emission band, which gets fainter and more narrow towards 
the centre, again due to a summation effect, which is stronger 
for the case of a line of sight tangential to the ring than perpen- 
dicular. 

At an inclination angle of 60 °, the same features remain, but 
the two semilunar areas of maximal emission get connected via 
thin and faint bows, due to the different line of sight to the ring- 
like structure. This effect further increases at / — 30 °. Now, the 
whole ring of the main dust grain contributors is visible. We see 
here the intersection of the sublimation sphere of these kind of 
grains with the torus. 

In the face-on case (z = °), we get the expected ring-like 
structure. In the central area no emission is present, because 
in our modelling the funnel of the dusty torus itself is free of 
any material. The maximum surface brightness does not result 
from the innermost part of the torus - although these dust grains 
possess the highest temperatures - because of the summation 
effect already mentioned above. Only very little dust can ex- 



ist here because of the density model itself and because the 
sublimation radii of the main dust contributors lie further out. 
Summing up along the line of sight results in the shown distri- 
bution of surface brightness. 

3.5. Wavelength study for surface brightness 
distributions 

In Fig. 10 we show surface brightness distributions at different 
wavelengths for the face-on and the edge-on view on our mean 
Seyfert torus model. Images for the same inclination angles are 
given with the same logarithmic colour scale and a dynamic 
range between the global maximum at the respective inclina- 
tion angle and 10~ 3 % of this value. Again, only re-emission 
spectra of the dust configuration are shown. The direct radia- 
tion of the energy source is neglected. When we would take this 
into account, the centre would always yield the highest bright- 
ness in the surface brightness distribution. The wavelengths are 
chosen as follows: We start with 7.22 (xm, located outside the 
silicate feature at 9.7 |im, but closest to 8 |xm, where the wave- 
length range of the MIDI (Leinert et al. 2003) instrument starts. 
The next is the wavelength of the feature itself, then between 
the two silicate features (at 13.18 |im) at the upper end of the 
MIDI wavelength range and finally we chose the overall maxi- 
mum of the SEDs (30.2 ^im). 

At wavelengths smaller than approximately 2 |.im, re- 
emission of dust can be neglected and pure emission of the cen- 
tral source remains, attenuated by the dust. The surface bright- 
ness distributions at varying wavelengths highlight explicitly 
the different temperature domains. According to Wien's law, 
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Fig. 10. Surface brightness distributions of our standard model for different wavelengths from left to right 
(7.22 urn, 9.7 |im, 13.18 (Am, 30.2 urn). Maps are given in logarithmic scale for inclination angles of ; = 0° (upper row) and i = 90° 
(lower row). See text for a detailed description. 



the maximum surface brightness at a wavelength of about 7 (Am 
is emitted by dust of a temperature of around 400 K. This area 
is concentrated to the innermost part of the dust distribution, 
caused by the very steep radial temperature curve there. Dust 
heated to about 100 K - which is reached in the flattening part 
of the temperature distribution - has its maximum radiation at a 
wavelength of around 30 (Am. Therefore, emission can be seen 
from a large fraction of the torus volume. The same behaviour 
is also visible in the case of an inclination of 90°. While for 
small wavelengths basically the x-shaped structure is visible, 
at longer wavelengths one can see the whole torus. This is also 
due to the fact that the extinction curve drops off very steeply 
and, therefore, we expect emission from all over the torus body. 



3.6. Implementation of a radiation characteristic 

A further step towards a more physical model was the intro- 
duction of a radiation characteristic of the accretion disc (see 
also Manske et al. (1998)). We wanted to check if this yields 
less steep radial temperature distributions, helping us to ex- 
plain the increasing depth of the silicate feature seen in corre- 
lated flux measurements performed by MIDI (Jaffe et al. 2004). 
In an extended optically thick accretion disc, a significant part 
of the radiation arises from the surface. Therefore, in a simple 
approach, we implemented a radiation characteristic following 
the same | cos(0)| law for all wavelengths. After emission of a 
photon package from the central source, the direction of emis- 
sion is chosen by chance, but in accordance with this radiation 
characteristic. 
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Fig. 11. Comparison of the temperature distribution of our standard 
model (solid lines) with the model with an anisotropically radiating 
accretion disc (dotted lines). The respective pairs of radial tempera- 
ture distributions for the smallest silicate grain component are shown. 
From left to right: within the equatorial plane, for an inclination angle 
of 30° and an inclination angle of 20°. 



In Fig. 11, the differences in the temperature distributions 
of the smallest silicate dust component are visualised. Here, 
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Fig. 12. SEDs at inclination angles of ° and 90 ° for our standard model with an isotropically radiating accretion disc (given by the solid line) 
and for the same model but with an accretion disc with a \cos(6)\ - radiation characteristic, without zoom (dotted line) and zooming into the 
inner 10 pc in radius (dashed line). 



one can directly see that compared to the isotropic radiating 
case, the | cos(#)l characteristic for the photon emission leads to 
hotter areas near the funnel walls, while the temperature within 
the equatorial plane gets lower. This is understandable, because 
the midplane is mainly heated by indirect reradiation of the dust 
in this case. 

In Fig. 12, SEDs of our standard model with an isotropi- 
cally radiating accretion disc (given by the solid line) are com- 
pared to the case of an accretion disc with a |cos 6\ - radiation 
characteristic (dotted line). The dashed line corresponds to the 
model with a radiation characteristic, but zoomed into the cen- 
tral lOpc in radius and will be discussed in Sect. 3.10. As the 
maximum dust temperature is smaller when implementing a ra- 
diation characteristic, the left branch of the IR bump is slightly 
shifted to longer wavelengths. The dotted line coincides here 
with the dashed line. With increasing inclination, this effect 
gets less important due to increasing extinction along the line 
of sight. In disagreement with Manske et al. (1998), we cannot 
find any change of the depth of the silicate feature, when com- 
paring the SED with an isotropically radiating accretion disc 
with the anisotropic case. This is caused by the different geom- 
etry in the innermost part of the torus and other values of the 
optical depth of our modelling compared to their approach. The 
cusp of the torus, which is only present in our modelling leads 
to a more pronounced silicate feature in emission for small in- 
clination angles. 

In the second row of Fig. 9, surface brightness distributions at 
different inclination angles are shown for the case of our stan- 
dard torus model, illuminated by the accretion disc described 
above. Here, one can directly see that the radiation character- 
istic leads to an enhanced heating of the funnel walls. This is 



best seen when comparing the images at an inclination angle 
of °. The area of maximum emission increases for the case of 
an implemented radiation characteristic, because more energy 
is emitted into the relevant solid angle and can, therefore, heat 
the funnel efficiently up to a larger height. In contrast to this, the 
equatorial plane can only be heated by reemitted photons and, 
therefore, only smaller temperatures are reached. This effect is 
visible directly when looking at the image for the edge-on view 
to the torus. The midplane appears as a narrow absorption band. 

3.7. Dust mass variation study for SEDs 

To study the impact of changing the dust mass, which is 
enclosed in the torus, we started with our standard model 
(given by the solid line in Fig. 13) with an optical depth of 
i"9.7(»m = 2.0, doubled the mass (dotted line, T9. 7kim =4.0) and 
halved it twice (dashed line with T9.7 ^ m = 1 .0 and dashed-dotted 
line with T9. 7[lm =0.5). The resulting SEDs are shown in Fig. 13 
for the SeyfertI and Seyfertll case with inclination angles of 
and 90 respectively. As the optical depth is given by the 
integral of the product of the dust density and the mass extinc- 
tion coefficient along the line of sight, it scales linearly with the 
enclosed dust mass: 

Text = J PAKzxlds. (10) 

i 

When we choose a line of sight within the dust-free hole of the 
torus - corresponding to the Seyfert I case - the whole temper- 
ature range of the dust in the torus can be seen. With increasing 
optical depth, more and more photons are absorbed within a 
decreasing volume near the funnel of the torus and reemitted in 
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Fig. 13. SEDs at inclination angles of 0° and 90° for different optical depths along an equatorial line of sight: Tc,.7 (lm = 4.0 (dotted line), 
Tq.i llm = 2.0 (solid line - our standard model ), Tc, 7 klra = 1.0 (dashed line) and T 9 7klm = 0.5 (dashed-dotted line). 



the infrared wavelength range. This leads to a rise of the tem- 
perature of dust near the funnel. But due to the high optical 
depths, less photons reach the outer part of the torus and, there- 
fore, the temperature drops. Altogether, a steeper temperature 
distribution in radial direction results. Smaller minimum tem- 
peratures lead to a shift of the Rayleigh-Jeans branch to longer 
wavelengths and larger maximum temperatures lead to a shift 
of the Wien branch to shorter wavelengths. Altogether, the IR 
bump gets broader. The higher temperatures in the inner re- 
gion also cause more flux at short wavelengths, while the larger 
amount of cold dust produces more flux at longer wavelengths. 
The global maximum of the IR bump is shifted to higher wave- 
lengths because of the smaller temperatures at the point of the 
maximum of the dust density. 

For an inclination angle of 90°, which corresponds to a 
Seyfertll-like case, there is always dust on the line of sight. 
Only the outer and, therefore, the coldest part of the torus can 
be seen directly. That is why we can use the reasoning for 
the SeyfertI case for the long wavelength region. The inner- 
most (and therefore hottest) part of the torus can only be seen 
through an increasing amount of dust. The higher the extinc- 
tion the less we can look into the torus and the less hot ar- 
eas of the torus can be seen. Therefore, the left branch of the 
IR bump is shifted to higher wavelengths when the extinction 
increases. The 9.7 \im silicate feature changes from emission 
over partial absorption to total absorption, comparable to large 
inclination angles in the study shown in Sect. 3.3. As discussed 
above, increasing dust mass affects the optical depth linearly. 
When t = 1 is reached, the torus gets optically thick, and, there- 
fore, the feature changes from emission to absorption. Between 
these two cases, a smooth transition region is visible, where the 
feature can be seen in partial absorption and emission, called 



self-absorption (Henning 1983). Intersecting the r = 1 line with 
the optical depth plotted against the wavelength though dis- 
plays the expected behavior of the silicate feature in the spec- 
trum. The same holds for the second silicate feature at 18.5 (im 
but in alleviated form as it is much weaker and partially cov- 
ered by the global maximum of the IR bump. 
The study shows that the SEDs are very sensitive to changes 
in the dust mass. Already a change by a factor of four changes 
the feature from emission to complete absorption (see solid line 
compared to dashed-dotted line in Fig. 13). Therefore, the sil- 
icate feature is a quite sensitive tool to distinguish the optical 
depth along the line of sight. 

3.8. Dust mass variation study for surface brightness 
distributions 

To study the impacts of changing the optical depth along the 
line of sight on surface brightness distributions, we used the 
same dust masses as in Sect. 3.7 for the case of spectral en- 
ergy distributions to gain comparability, but doubled the high- 
est mass again. This means that coming from our standard 
model, we doubled the enclosed dust mass twice and halved 
it twice. Again, we chose the same wavelength of 13.18 (im. 
The results are shown in Fig. 14. The different optical depths 
(T9.7 Mm = °- 5 > T 9.7 = 1 -0, t 9 . 7 m = 4.0, t 9j ^ m = 8.0) are shown 
in different columns, while the inclination angle changes are 
shown in different rows, starting from 0° in the upper row 
to 90 ° in the lower row. Our standard model is omitted here, 
please compare to Fig. 9. All of the images are given in a lin- 
ear colour scale, with a dynamic range between zero and the 
maximum value of the respective map in order to show the 
changes in the distribution of maximum surface brightness for 
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Fig. 14. Surface brightness distributions at a wavelength of A = 13.18 um for models with different optical depths. The values for the optical 
depth are 79.7^ = 0.5 (left column), 79.7^ = 1.0, Tgj flm = 4.0 and T 9 .7^ m = 8.0 (right column) and inclination angle changes from (upper 
row) to 90° (lower row) in steps of 30°. Each step in r means doubling the included dust mass. The images with 79.7^ = 2.0 are omitted 
(please compare to Fig. 9). All images are given with a linear colour scale reaching from zero to the maximum of the respective image. 



the different optical depths. The absolute values for the pixels 
are given in Fig. 15, which shows a cut through the maps of 
the dust mass variation study from the centre along the positive 
x-axis for an inclination of 0° and 90°. The fluxes are given 
in Jansky per square milliarcsecond for an assumed distance of 



the object of 45 Mpc. For the case of °, all graphs show the 
same trend: coming from the centre, one can see a sharp rise 
of surface brightness with a sharp peak and a flatter decrease in 
the outer part. Between 4 and 5 pc, all of the curves meet into 
a single curve, independent of the optical depth. The latter is 
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i=0° i=90° 




Fig. 15. Cut through the surface brightness distributions fi v at a wavelength of A = 13.18 |im along the positive x-axis for inclination angles of 
0° and 90°. The thick solid line corresponds to the case with T 97(lm =0.5, the thin solid line to T9 7(ml = 1.0, the dotted line to our standard 
model with T9. 7(lm = 2.0, the thick dashed line to 79.7,1m =4.0 and the thin dashed line to T9. 7(ml = 8.0. The fluxes shown are given in Jansky per 
square milliarcsecond for an assumed distance of the object of 45 Mpc. 



due to the fact that our model features very steep funnel walls. 
As indicated by the temperature distribution, we expect most 
of the emission for this wavelength next to the torus boundary. 
The intersection of all curves then indicates the distance from 
the centre, where the density has decreased that much that a 
doubling of it does not change the emitted flux anymore and 
the expected surface brightness distribution is independent of 
the optical depth within our sample. By increasing the optical 
depth, the maximum is shifted towards the centre and to higher 
values and the curves get steeper. As already discussed in the 
previous section, the higher the optical depth is the steeper the 
drop of the temperature distribution in radial direction is, be- 
cause the photon mean free path length decreases and most of 
the photon packages get absorbed within a reduced volume of 
the dust torus (steepening of the inner part) and, therefore, only 
reemitted photons in the IR reach the outer region (flattening 
of the outer part). To estimate the surface brightness, one has 
to take into account all emission by dust until t = 1 is reached. 
This can also be directly seen in Fig. 14 in the decreasing size 
of the region of maximum emission. The apparent shift of the 
inner radius of the torus is caused by the increase of the atmo- 
sphere around the dense core of the torus or the growth of the 
denser core, respectively. 

At the other inclination angles one can again see that 
the area of maximal emission changes. At smaller optical 
depths, it is a diffuse, wide distribution within the torus vol- 
ume. Increasing the optical depth leads to a spreading, which 
is more located in the direct vicinity of the inner funnel bor- 
der. This leads to the sharply emerging inner borders of the 
funnel, which are now visible in the typical x-shaped form. 



These appearances have been seen in a similar dust config- 
uration around a double system, the so-called Red Rectangle 
(Tuthill et al. 2002), which can be explained by a dusty torus 
as well (Men'shchikov et al. 2002). It is located in a distance of 
only 330 pc, much closer than Seyfert galaxies and can there- 
fore be resolved. 



At higher inclination angles, two effects have to be taken 
into account. With increasing amount of dust along the line of 
sight, first of all, we expect a larger surface brightness, but the 
extinction along the line of sight increases as well. The value 
of the maximum of the cut along the positive x-axis of the 
maps with the highest optical depth is not the brightest any- 
more. It decreases with increasing inclinations. In Fig. 15b the 
case for an inclination angle of 90 is shown. Here, the model 
with the highest dust mass (thin dashed line) has the smallest 
surface brightness and the curve of the second heaviest torus 
(thick dashed line) has also shifted significantly to smaller val- 
ues. The shape of these graphs also changes when changing 
from the Seyfert I (i = °) case to the Seyfert II (i = 90 °) case. 
As the directly illuminated parts of the funnel are not directly 
visible anymore, the decreasing curve further away from the 
centre decreases faster. With increasing optical depth, the outer 
part of the torus gets colder and colder, which can be seen by 
the decreasing flux and the steepening of the curves there. In 
the innermost part, the curves flatten, which is due to the emis- 
sion band of the innermost ring of the funnel. As less photons 
reach the outer part, statistic gets poorer there and therefore the 
graphs get more noisy. 
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Fig. 16. Optical depth within the equatorial plane for different dust models: a Increasing the exponent of the number density distribution of 
grain sizes, from -3.5 (our standard model, given by the solid line) over -3.0 (dotted line) and -2.5 (dashed line) to -2.0 (dashed-dotted line), 
b Changing the range of grain sizes, from 0.005 . . . 0.25 |im (our standard model, given by the solid line) to 0.005 ... 10 urn (dotted line) and to 
0.001... 10 um (dashed line). 



3.9. Variation of dust properties 

As already pointed out in Sect. 2.3, the dust composition in 
the extreme environment of an AGN may be altered compared 
to interstellar dust in our own galaxy. The main effects here 
are sputtering due to interaction of the dust grains with the hot 
surrounding gas, the presence of the strong radiation field of 
the central source, destruction of dust grains due to shocks, etc. 

All of the mentioned effects mainly impact the population 
of small grains. Therefore, some authors (e. g. Laor & Draine 
1993; Maiolino et al. 2001) claim that within an AGN dust 
model, small grains should be depleted. In the following sec- 
tion, we test this by first of all changing the exponent of the 
number density distribution of grain sizes and later by widen- 
ing the grain size distribution in both directions. By doing this, 
a transition is made between the standard MRN galactic dust 
model - which we have adopted as our standard dust model 
- and the dust model proposed by van Bemmel & Dullemond 
(2003). In this section, we are especially interested in the be- 
haviour of the 9.7mn silicate feature, as it was claimed that the 
depletion of smaller grains may solve the problem that the fea- 
ture is never seen in emission in AGN. On the other hand one 
needs to check whether other constraints are violated by this 
change of the dust model, e. g. the appearance of the feature in 
absorption, as observed. 

3.9.1 . Changing the exponent of the number 
distribution of grain sizes 

Different exponents of the number density distribution of the 
grains were tested by van Bemmel & Dullemond (2003) for 
their torus model. They finally adopt a somehow more shallow 
distribution with an exponent of -2 to be their standard dust 
model. 

In this section, we also check different slopes of the dis- 
tribution for our standard torus geometry (see Sect. 2.2). We 
start with our standard dust model with a slope of -3.5 (given 
by the solid line in Fig. 16a and Fig. 17) and test shallower 



distributions by changing the exponent in steps of 0.5 up to 
-2.0 (dashed-dotted line) used by van Bemmel & Dullemond 
(2003). Due to renormalisation of the dust density after chang- 
ing the exponent, this means that a flatter distribution leads to 
more large grains compared to small grains. This produces a 
flatter course of the optical depth (see Fig. 16a) at small wave- 
lengths - compared to our standard model - because of the less 
efficient heating of larger grains, which dominate the dust com- 
position more and more during the variation. Consequently, the 
temperature distributions flatten. 

In the SEDs (see Fig. 17a), this can be seen by the de- 
creasing fluxes at the rise of the IR-bump at small wavelengths 
and the increasing fluxes at wavelengths around the maximum 
(20 - 30p,m) of the infrared bump (see Fig. 17a). The strati- 
fication of dust grains concerning sizes and composition and 
the increasing extinction caused by the large grain population 
(with smaller sublimation radii) lead to larger sublimation radii 
for the other components. In contradiction to the modelling of 
van Bemmel & Dullemond (2003), a significant change of the 
silicate feature cannot be found in our simulations. For the case 
of an inclination angle of 90° (see Fig. 17b), any change is sim- 
ply due to the fact that there is too little dust left that has tem- 
peratures higher than about 300K - the temperature which is 
required for the emission feature at 9.7mn. 

3.9.2. Changing the width of the grain size distribution 

In a further dust parameter study, we test the effects of a broad- 
ening of the grain size distribution. The corresponding courses 
of the optical depths with wavelength within the equatorial 
plane are plotted in Fig. 16b and the resulting SEDs in Fig. 18b. 
Starting from the MRN model (given by the solid line), we first 
extend the grain size range towards larger grains (dotted line) 
up to lOum size and later additionally towards smaller grains 
down to 0.001 Jim (dashed line). By doing this, we finally get 
the dust model used by van Bemmel & Dullemond (2003). 

As can be seen from Fig. 16b, only minor changes of the 
optical depth are visible. Conspicuous is the smaller relative 
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Fig. 17. Consequences on the SEDs at inclination angles of 0° and 90° of changing the exponent of the dust grain number density distribution. 
The exponent changes from -3.5 (solid line) in steps of 0.5 to -2.0 (dashed-dotted line). 
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Fig. 18. Consequences on the SEDs at inclination angles of 0° and 90° of changing the width of the dust grain size distribution from 
0.005 . . . 0.25 urn (solid line) over 0.005 . . . 10 (im (dotted line) and finally to 0.001 . . . 10 urn (dashed line). 



depth of the silicate feature towards shorter wavelengths after 
including also large grains (dotted line). This can also be seen 
in the SEDs, displayed in Fig. 18b. While the differences for an 
inclination angle of 0° are nearly negligible, the reduced rela- 
tive depth of the silicate feature towards smaller wavelengths is 
clearly visible in the SED for edge-on view (z = 90°). In the last 
step of our study - the inclusion of pronounced slightly smaller 



grains - the effects of extending the size distribution towards 
larger grains is partly compensated. 
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Fig. 19. Dependence of SEDs on the beamsize of the observing device. The solid line corresponds to our standard model, the dotted line to the 
resolution obtainable with a single-dish telescope (100 mas) for a distance to the object of 45 Mpc and the dashed line to the resolution of the 
MIDI interferometer (10 mas). The left panel shows SEDs for an inclination angle i = ° and the right panel for i = 90°. 



3.10. Zooming into the torus 

Another simulation series we made was to investigate the im- 
pact of zooming into the torus resembling a certain aperture of 
the telescope. To do this, we again started from our standard 
model and zoomed into the central lOpc and lpc in radius. 
This corresponds to a resolution of the observing instrument 
of 100 mas and 10 mas for an object at a distance of 45 Mpc. 
These are typical values for the largest single-dish telescopes 
(100 mas) and the first interferometers in this wavelength range 
like MIDI (10 mas). The results of this study are shown in 
Fig. 19. In terms of surface brightness distributions compare 
to Fig. 9. 

When we decrease the aperture, we concentrate more 
and more on the central parts of the torus emission. As our 
temperature distribution is very steep, it is expected that cold 
dust in the outer part of the torus is excluded from the field 
of view. Therefore, the smaller the aperture, the less total flux 
we obtain, which is taken from the large wavelengths part of 
the spectrum. For the case of 10 mas resolution, only the inner 
pc in radius is visible. This is even less than the sublimation 
radius of the smallest silicate grains (approximately 1.3 pc). 

For the silicate feature, one notices a decrease of the rela- 
tive depth rather than the expected increase, as seen with MIDI 
(compare to Fig. 2 in Jaffe et al. (2004)). This might be due 
to a too steep radial temperature distribution of our models. 
Introducing a cloudy structure, which leaves unobscured chan- 
nels up to the outer part of the torus, may lead to a direct heat- 
ing of clouds further outside and therefore to a broader tem- 
perature distribution. Zooming into the centre then means for 



the Seyfert II case to exclude areas, where emission features are 
produced and a deepening of the silicate absorption feature is 
expected to take place. 

4. Comparison with observations 

In this section, simulated spectral energy distributions are com- 
pared to available observational data for the case of large aper- 
ture SEDs and high spatial resolution SEDs. 

4.1. Comparison with large aperture spectra of type I 
galaxies 

In Fig. 20 we compare the resulting SEDs of our standard 
model (different line styles) for an axial view on the torus 
(typically, Seyfert I galaxies have inclination angles between 
Wand 20°) with varying enclosed dust masses (as described 
in Sect. 3.7 and shown in Fig. 13) to different sets of mean large 
aperture data for type I objects. The plus signs refer to a mean 
spectrum, extrapolated from 29 radio quiet quasars, provided 
by Elvis et al. (1994), the squares correspond to 16 Seyfert 
galaxies, reported by Granato & Danese (1994). Triangles and 
diamonds represent radioquiet type I PG quasars with bolo- 
metric luminosity larger and smaller than 10 12 Lq respectively 
(Sanders et al. 1989). All of these data sets are scaled in order 
to get fluxes comparable to our simulations. 

As can be seen from this, our models compare well to the 
overall shape of the UV to FIR wavelength range of these mean 
AGN spectra. An exeption is the pronounced feature at 9.7 urn. 
As already discussed before, it has never been observed in 
emission in SEDs of Seyfert galaxies - as well not in SEDs of 
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Fig. 20. Comparison of our standard model (different line styles cor- 
respond to varying dust masses) with large aperture spectra for type I 
objects. Plus signs refer to 29 radio quiet quasars (Elvis et al. 1994), 
squares to 16 SeyfertI galaxies (Granato & Danese 1994) and dia- 
monds and triangles to PG Quasars (Sanders et al. 1989) with lumi- 
nosity smaller (diamonds) and larger (triang les) than 10 12 L o . For a 
typical SeyfertI galaxy (10° < i < 20°), the silicate feature looks 
much less pronounced (see Fig. 8a). 



4.2.1. NGC1068 

NGC 1068 is a nearby Seyfertll galaxy at a distance of 
14.4 Mpc (Galliano et al. 2003). Therefore, it has been stud- 
ied extensively and we were able to find at least some of the 
needed simulation parameters in literature. The procedure de- 
scribed above led to the parameters summarised in Table 2. The 
inclination angle of i — 85 ° was chosen in agreement with the 
tilt of the outflow cone of ionised gas reported by Crenshaw 
& Kraemer (2000). For our comparison, we use a compilation 
of high resolution data done by M. A. Prieto, given in Table 3. 
Fig. 21 shows the result of the adaptation of the model given 
in Table 2 to this data set. The line corresponds to a zoomed-in 
view of the torus by a factor of 10, which resembles the aper- 
ture of 0.2"in diameter of most of the datapoints (given by the 
triangles) in Fig. 21. The diamonds correspond to a slightly 
larger aperture of 0.27". The L-band flux, which corresponds 
to an aperture of 0.7"in diameter is given by the square sym- 
bol. In addition, we added the total flux (dashed-dotted line) - 
corresponding to an aperture of 0.6 x 0.6" - and the correlated 
flux spectrum for a baseline of 78 m (dotted line) - correspond- 
ing to an aperture in diameter of approximately 30 mas (Jaffe 
et al. 2004). As can be seen in this figure, observations are in 
good agreement with our model except of the wavelength re- 
gion around the 9.7 \im silicate characteristic. As already men- 
tioned before, zooming into our model torus does not reproduce 
the trend of MIDI-observations, which means a deepening of 
the silicate feature. We therefore adapted the depth of our fea- 
ture to the depth observed in correlated flux measurements of 
MIDI with a baseline of 78 m. This was also necessary in order 
to describe the datapoints at the smaller wavelengths. 



higher spectral resolution. The introduction of a clumpy struc- 
ture of the dust distribution, which is also a more physical so- 
lution, might help to remove these large emission features of 
our models. A comparison with high spatial resolution SEDs 
of type II sources is given in the next section for some selected 
objects. 

4.2. Comparison with special Seyfertll galaxies 

To specify as many of our simulation parameters as possible for 
these special objects, we started with a literature study. Some 
of the remaining variables could be chosen according to con- 
ditional relations as mentioned in Sect. 2.2. The value of the 
enclosed dust mass was assumed such that the relative depth 
of the 10 ^m-feature appears comparable to the observations of 
the correlated flux with a baseline of 78 m done with MIDI for 
the case of NGC 1068 and comparable to the depth of the fea- 
ture in the MIR spectrum of Roche et al. ( 1 99 1 ) for the Circinus 
galaxy. As zooming into our torus does not reproduce the deep- 
ening of the silicate feature as seen by MIDI, this seems to be 
a reasonable procedure and explains the deviations concerning 
wavelengths around the feature for the case of NGC 1068. It 
results in an optical depth in the equatorial plane at 9.7 (im of 
t = 3.3 for NGC 1068 and t = 3.9 for the Circinus galaxy. The 
free parameters finally were chosen to obtain a reasonable dust 
distribution and a good comparison with the data. 



Table 2. Parameters used for the simulation of NGC 1068. 



Parameter 


Value 


Reference 


Mbh 


8.3 • 10" M Q 


Greenhill et al. (1996) 


R c 


24 pc 


Gallimore & Matthews (2003) 


M, 


6.5 • 10 8 Mq 


Thatteetal. (1997) 




5.42- 10 10 L o 




Alisc / ^Edd 


20% 




Rj 


3pc 




M dust 


7.8 • 10 4 Mq 




T 9.7 |im 


3.3 




^out 


72 pc 




v t 


164 km s- 1 




r 


0.5 




d 


14.4 Mpc 


Galliano et al. (2003) 


i 


85° 


Crenshaw & Kraemer (2000) 



4.2.2. Circinus 

The Circinus galaxy is one of the nearest spiral galaxies, which 
harbors an Active Galactic Nucleus. As it is the second ex- 
tragalactic object, which has been observed by MIDI and also 
with NAOS/CONICA (Prieto et al. 2004), we also carried out 
a comparison with one of our models. From the first results of 
these observations (Prieto et al. 2004) and also from previous 



Marc Schartmann et al.: Towards a 



i=85° 

10 36 r , , 1— 




A [/xm] 

Fig. 21. Comparison of one of our models (described by the parame- 
ters given in Table 2 and plotted with a solid line) to a set of data (given 
in Table 3) from high spatial resolution observations of NGC 1068, 
compiled by M. A. Prieto. The various symbols correspond to slightly 
different apertures. The dashed-dotted line gives the total flux spec- 
trum obtained with MIDI and the dotted line corresponds to the corre- 
lated flux spectrum obtained with MIDI for a baseline of 78 m (Jaffe 
et al. 2004). 

Table 3. High resolution data of NGC 1068. The compilation was 
done by M. A. Prieto. The references are R98 (Rouan et al. 1998), P04 
(Prieto, private communication), M03 (Marco & Brooks 2003) and 
BOO (Bock et al. 2000). 



Wavelength 


Band 


Flux 


Aperture (dia) 


Reference 


[im 




Jy 


arcsec 




1.30 


J 


5e-4 


0.20 


R98 


1.70 


H 


6e-3 


0.20 


R98 


2.20 


K 


5.6e-2 


0.27 


P04 


3.79 


L 


1.5 


0.70 


M03 


4.47 


M 


2.5 


0.27 


P04 


7.99 


N 


4.23 


0.20 


BOO 


9.99 


N 


3.33 


0.20 


BOO 


11.99 


N 


7.1 


0.20 


BOO 


24.98 




9.4 


0.20 


BOO 



modelling of the dusty torus by Ruiz et al. (2001), we know 
that it harbors a relatively small sized toroidal dust distribu- 
tion. Therefore, we used a scaled down version of our standard 
model described above. The simulation was done with the pa- 
rameters summarised in Table 4. 

Fig. 22 shows the result of the comparison of this model 
(given by the solid line plot) with data taken with the NACO 
camera at the VLT (Prieto et al. 2004). Various symbols rep- 
resent different apertures (see Table 5). Fluxes, corrected for 
foreground extinction by Ay = 6 mag, are used. 
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Fig. 22. Model comparison with high spatial resolution data of the 
Circinus galaxy taken with the NACO camera at the VLT. Data cour- 
tesy of Prieto et al. (2004). The data was corrected for foreground ex- 
tinction by Am = 6 mag. The various symbols correspond to different 
apertures (see Table 5). 



Table 4. Parameters used for the simulation of the Circinus galaxy. 



Parameter 


Value 


Reference 


Mbh 


1.7 • 10" M Q 


Greenhill et al. (2003) 


Rc 


10 pc 




M, 


2.0 • 10 8 Mq 




Aiisc 


8.4- 10 9 Lq, 




-^disc/^Edd 


15% 






3.5 pc 




M dust 


2.0 • 10 4 Mq 




Tg.7 urn 


3.9 




^out 


30 pc 




v t 


140kms-» 




r 


0.5 




d 


4.2 Mpc 


Freeman et al. (1977) 


i 


90° 


Crenshaw & Kraemer (2000) 



Table 5. High resolution data of the Circinus galaxy. The data was 
corrected for foreground extinction by A v = 6 mag. Data courtesy of 
Prieto et al. (2004). 



Wavelength 


Band 


Flux 


Aperture (dia) 


\im 




mJy 


arcsec 


1.30 


J 


S7.9 


0.38 


1.70 


H 


15.0 


0.10 


2.20 


K 


36.0 


0.38 


2.42 




54.0 


0.38 


3.79 


L 


480.0 


0.38 


4.47 


M 


2200.0 


0.38 


11.80 


N 


15000.0 


1.00 
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The graph corresponds to a torus model with an inner radius 
of approximately 0.9 pc and a luminosity of the central source 
of 15% of the Eddington luminosity. The modelling leads to 
a maximum temperature lower than the sublimation tempera- 
tures of the dust grains. Models with smaller inner radii and, 
therefore, higher maximum temperatures can be excluded. For 
models with higher maximum temperature, the rise of the IR 
bump at small wavelengths is moved towards smaller wave- 
lengths. In order to avoid this, we would have to increase the 
amount of dust, leading to a deeper silicate feature, in contra- 
diction to available observations. This shows that sublimation 
temperatures of the dust grains are only reached in some of the 
torus configurations. Our model fits well to the observations, 
except of the small wavelength part of the SED. This excess is 
not due to torus emission. The increased fluxes are most proba- 
ble caused by light from the central source, scattered above the 
torus by material (dust and electrons) within the opening of the 
torus. 

Our attempts to model observational data showed once 
more the inherent ambiguity: The same SED can be modelled 
with different parameter sets or even models. Therefore, addi- 
tional information from observations is needed in order to be 
able to distinguish between different torus models. 

5. Conclusions and outlook 

In this paper we present new radiative transfer simulations for 
dust tori of Active Galactic Nuclei. The dust density distribu- 
tion and the geometrical shape of the torus model result from a 
hydrostatic equilibrium in which gravitational and centrifugal 
forces of the nuclear stellar distribution and the central black 
hole are balanced by pressure forces due to the turbulent mo- 
tion of the clouds. Spectral energy distributions (SEDs) and 
surface brightness distributions are obtained by a fully three- 
dimensional treatment of the radiative transfer problem. 

After fixing as many parameters of the model as possible 
by observational constraints, we undertake an extensive param- 
eter study by varying dust properties, dust masses and other 
factors. It shows that our modelling is able to explain mean 
observed spectral energy distributions (SEDs) for classes of 
Active Galactic Nuclei as well as for individual objects. 

Obviously, the innermost part of the tori close to the cen- 
tral energy source is crucial for the determination of SEDs in 
the mid-infrared wavelength range. Specifically, we found that 
any realistic model has to incorporate the fact that sublima- 
tion radii for varying dust grain sizes and dust composition dif- 
fer considerably. Using mean dust characteristics according to 
the MRN dust model (see Sect. 2.3) instead leads to incorrect 
temperature distributions and SEDs. Even with this refinement, 
the model cannot explain the recent interferometric observa- 
tions of the MIDI instrument, as reported by Jaffe et al. (2004), 
which show a deepening of the silicate absorption feature in 
the correlated flux - effectively sampling the central few par- 
sec of the dust distribution. In order to reproduce these results, 
a temperature distribution shallower than in our model would 
be required. We find that introducing a | cos 6\ radiation char- 
acteristic for the primary radiation source does not yield sig- 
nificant changes. In a next step, we will introduce a clumpy 



structure of our dust configuration. Gaps within the dust distri- 
bution should then lead to direct heating of clumps in the outer 
part of the torus and thus producing a shallower mean tempera- 
ture distribution. Biasing the grain size distribution of the dust 
towards larger grains as proposed by several authors (e.g. Laor 
& Draine 1993; Maiolino et al. 2001) does not yield a reduc- 
tion of the emission feature of type I objects for the case of our 
dust density distribution. 

With the present examination of the Turbulent Torus Model 
for dust distributions in AGNs, we made a first step towards 
the introduction of more physics into the modelling of the torus 
emission. We demonstrated that a model based on physical as- 
sumptions reproduces the SEDs both on kiloparsec and par- 
sec scales as well as previous models in which the shape and 
size of the dust distribution are free parameters {e.g. Pier & 
Krolik (1993); Granato & Danese (1994); Manske et al. (1998); 
Nenkova et al. (2002); van Bemmel & Dullemond (2003) ). 

As pointed out in Sect. 2.1, this kind of stationary mod- 
elling cannot account for all physical effects present in such 
an environment. In a further step, hydrodynamic simulations in 
combination with radiative transfer calculations should enable 
us to take at least some of these uncertainties into account. The 
main effects are continuous energy feedback by supernovae ex- 
plosions and mass injection by stellar winds. 

With the first interferometric observations at hand, how- 
ever, the problem becomes more challenging: the wide range 
in observable scales from hundreds of parsec down to one par- 
sec requires not only to reproduce global SEDs of the torus 
but also to describe its shape, structure and temperature strat- 
ification correctly. In addition, these detailed observations are 
an essential step to resolve the ambiguities that are present in 
current global models of dusty tori in Active Galactic Nuclei. 
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